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This  report  presents  investigations  and  results  of  the  boundary-  and 
finite-element  techniques  for  solving  soil-structure  interaction  problems. 
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Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI 
(metric)  units  as  follows: 
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_ Multiply _  By _ 

feet  0.3048 

kips  (force)  4.448222 

kip  (force)-feet  1355.818 

kips  (force)  per  square  inch  6.894757 

pounds  (force)  per  square  inch  6.894757 


To  Obtain 

metres 
kilonewtons 
newton-metres 
megapascals 
kilopascals 


THE  APPLICATION  OF  BOUNDARY- ELEMENT  TECHNIQUES  FOR 
SOME  SOIL-STRUCTURE  INTERACTION  PROBLEMS 


PART  I:  PROBLEMS  AND  SOLUTIONS 

Introduction 

1.  Soil-structure  interaction  problems  result  when  the  behavior  of  the 
structure  and  the  surrounding  soil  are  interdependent,  and  the  solution  re¬ 
quires  an  analysis  of  both  the  structure  and  the  soil  in  a  compatible  manner. 
While  structures  are  often  satisfactorily  modeled  as  linearly  elastic,  homo¬ 
geneous,  and  isotropic  materials,  the  modeling  of  soils  is  extremely  complex. 
To  model  the  in-situ  behavior  of  soils,  one  must  make  gross  approximations  by 
using  experience  and  judgment,  and  these  evaluations  are  chiefly  based  on  the 
relative  importance  of  the  project  and  the  desired  accuracy.  The  complexities 
in  the  constitutive  relations  of  the  soil  continuum  are  enhanced  by  the  fact 
that  the  soil  has  been  deposited  in  nature  in  a  layered,  heterogeneous  manner. 
Since,  in  most  instances,  engineers  are  interested  only  in  the  behavior  of 
structures,  they  have  assumed  very  simplified  properties  of  the  soil  in  their 
design  considerations.  One  of  the  widely  used  models  for  soil-structure  in¬ 
teraction  problems  is  the  Winkler  spring  model  (Scott  1981)  for  analysis  of 
beams  on  elastic  foundations,  mat  foundations,  pavements,  pile  foundations, 
etc . 

Winkler  Model 

2.  Winkler,  in  1867,  proposed  that  the  deflection  of  the  soil  surface 
can  be  modeled  by  a  simple  equation,* 

p  --  kw  (1) 

where 

p  =  the  pressure  acting  on  the  soil  surface 

k  =  the  proportionality  constant,  known  as  the  subgrade  modulus  or  the 
modulus  of  subgrade  reaction 

w  =  the  deflection  of  the  loaded  region  on  the  surface 


The  unit  of  k  is  in  pounds  per  cubic  inch.  Normally,  simple  plate  bearing 
tests  are  conducted  to  determine  the  value  of  k  .  This  concept  has  been 
widely  used  by  engineers.  For  a  given  value  of  k  ,  Hetenyi  (1946)  solved 
many  problems  of  beams  on  elastic  foundation,  and  Westergaard  (1926)  solved 
problems  of  slabs  on  elastic  foundations.  When  the  question  of  the  value  of 
k  of  soil  was  raised,  Terzaghi  (1955)  offered  general  guidelines.  Matlock 
and  Reese  (I960)  have  widely  used  this  technique  for  solving  problems  of  lat¬ 
erally  loaded  pile  foundations.  Terzaghi  showed  that  though  linearly  elas¬ 
tic,  isotropic,  and  homogeneous  properties  for  soil  were  used,  the  modulus 
of  subgrade  reaction  depended  very  much  on  the  size  of  the  loaded  area. 

Vesic  (1961)  showed  that  the  value  of  k  is  influenced  by  the  stiffness  of 
the  beam.  Based  on  experiments  and  theoretical  concepts,  many  empirical 
formulas  arose  for  values  of  k  .  These  are  always  questioned  by  structural 
engineers  who  generally  demand  a  soil  value  of  k  for  their  input  for  the 
soil-structure-interaction  analysis.  The  absence  of  a  unique  value  for  k 
for  soil  can  be  easily  realized  from  the  following  examples,  even  though  one 
assumes  idealized  properties  such  as  linear,  elastic,  etc.,  for  the  soil 
cont inuum. 


Example  Problem 

3.  Consider  a  plane-strain  problem  with  a  strip  of  uniformly  loaded 
beam  resting  on  a  semi- inf in ite  soil  continuum.  If  the  beam  is  very  rigid  as 
compared  to  the  soil,  the  deflections  along  the  beam  are  uniform  while  the 
pressure  distribution  varies  from  infinity  to  a  finite  value  at  the  center. 
This  is  shown  in  Figure  1.  On  the  other  hand,  if  the  beam  has  very  low  rigid¬ 
ity  its  deflections  vary  from  a  maximum  value  at  the  center  to  a  smaller  value 
at  the  ends,  as  shown  in  Figure  2.  In  each  case,  the  ratio  of  pressure  versus 
deflection  is  not  a  constant  at  the  soi l -structure  interface,  hence  the  non¬ 
uniqueness  of  the  value  of  k  is  demonstrated  here.  Realizing  these  prob¬ 
lems,  Pasternak  (1954)  developed  a  two-parameter  model  to  take  into  consider¬ 
ation  the  end  effects.  The  evaluation  of  these  parameters,  however,  again 
becomes  a  major  problem  confronting  the  geotechnical  engineer.  In  addition, 
it  can  be  shown  that  the  value  of  k  depends  on  the  depth  of  the  soil  contin¬ 
uum;  in  other  words,  it  is  influenced  by  the  existence  of  hard-rock  stratum 
and  the  depth  at  which  it  occurs. 
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Constant  Displacement 


Figure  1.  Pressure  distribution  on  the  interface 
of  a  rigid  beam 


Figure  2.  Deflection  on  the  interface  of  a 
flexible  beam 

4.  In  spite  of  these  limitations,  the  concept  of  the  constant,  or  even 
nonlinear,  k  is  used  by  engineers  for  analysis  of  moderately  simple  struc¬ 
tures.  However,  when  these  concepts  are  applied  to  analysis  of  major  struc¬ 
tures,  such  as  hydraulic  U-iock  structures,  a  better  md  more  accurate  deter¬ 
mination  of  the  soil  stiffness  based  on  its  const  i  tut  i  ve  t  •*;  it  ions  is 
necessary . 


Alternate  Methods 

5.  An  alternate  and  better  method  of  analysis  is  to  mod-  1  the  involved 
soil  continuum  in  its  entirety.  Closed  form  mathematical  solutions  to 


this  become  too  tedious,  even  for  relatively  simple  problems.  However,  elec¬ 
tronic  computers  made  it  possible  to  analyze  these  problems  numerically,  by 
using  methods  such  as  finite  difference,  finite  element,  etc.  Researchers, 
Duncan  and  Clough  (1971)  and  Vallabhan  and  Jain  (1972)  have  used  finite- 
element  methods  (FEM)  for  solving  soil-structure  interaction  problems.  These 
methods  become  cumbersome  as  the  number  of  unknowns  increases  rapidly,  and 
more  so  when  the  discretizations  of  the  continuum  are  altered  for  more  accu¬ 
rate  representations  of  the  soil  continuum.  Normal  discretizations  of  the 
soil  and  the  structure  result  in  a  fairly  coarse  discretization  for  the  struc¬ 
ture.  It  is  in  this  context  that  the  boundary-element  method  (BEM)  is  found 
advantageous  in  the  analysis  of  soil-structure  interaction  problems. 

6.  Two-dimensional  continuum  problems  are  represented  at  the  boundary 
as  line  elements  and  three-dimensional  problems  are  modeled  by  the  boundary 
surface  as  surface  elements,  i.e.,  BEM,  offer  a  reduction  in  the  dimension¬ 
ality  of  the  problem.  It  is  found  to  be  an  elegant  procedure  for  handling 
large  soil  continuum  with  elements  required  only  on  the  boundaries,  hence, 
attractive  to  the  analyst. 

Capabilities  of  the  Boundary-Element  Method 

7.  Engineers  are  seeking  parameters  to  represent  the  soil  stiffness 
underneath  the  structure,  with  simple  and  reasonable  assumptions  on  the  soil 
properties  such  as  Young's  modulus  E  and  Poisson's  ratio  v  of  the  soil 
medium.  Using  the  BEM,  the  stiffness  parameters  can  be  developed  in  a  matrix 
form,  without  depending  on  the  stiffness  of  the  structure,  size  of  the  soil- 
structure  interface  area,  and  the  distribution  of  the  loading.  The  emphasis 
of  this  research  report  is  to  investigate  the  boundary-element  techniques  to 
represent  a  soil  medium  in  such  a  manner  that  the  engineer  can  more  accurately 
use  it  in  his  soil-structure  interaction  problems. 


Mathematical  Aspects 

8.  The  mathematical  aspects  of  this  technique  are  given  in  detail  in 
Part  II.  Presented  here  is  an  overall  view  of  the  method,  highlighting  its 
advantages  over  domain-type  solutions.  The  dimensionality  of  the  problem  is 
reduced  by  one,  because  the  boundary  of  the  domain  alone  is  discretized  for 


analysis.  Therefore,  the  input  data  is  highly  simplified,  reducing  man  hours 
in  preparation  of  data.  Modeling  of  domains  extending  to  large  distances  are 
carried  out  very  efficiently.  Moreover,  this  method,  when  applied  to  a 
stress-concentration  problem,  yields  accurate  results  when  compared  to  other 
numerical  methods.  The  only  disadvantage  is  that  the  system  matrix  is  fully 
populated  and  numerical  efficiency  is  not  great  for  slender  regions,  unless 
the  domain  is  divided  into  regions  to  get  a  banded  type  matrix.  Overall,  the 
accuracy  and  efficiency  of  the  model  is  much  higher  than  those  of  the  other 
prevailing  methods  for  soil  structure  interaction  problems. 

Various  Methods  Used 

9.  In  this  report,  the  FEM  is  used  for  representing  the  U-lock 
structure,  and  the  BEM  is  used  for  modeling  the  soil.  Rectangular  plane- 
strain  elements  are  used  for  the  structure  in  the  FEM.  The  boundary  of  the 
soil  medium  is  discretized  as  constant  boundary  elements.  Using  a  static- 
condensation  procedure,  the  boundary-element  system  equations  are  reduced  to 
the  unknown  displacements  on  the  soil-structure  interface.  The  stiffness  of 
the  soil  system  is  thus  developed  on  the  interface  and  added  to  the  finite- 
element  structure  stiffness  matrix.  The  resulting  global  stiffness  is  solved 
for  stresses  and  displacements  of  the  structure,  including  the  interface.  The 
pressure  distribution  on  the  interface  is  then  calculated  from  the  boundary- 
element  equations,  using  the  interface  displacements.  Part  III  discusses  the 
theory  of  the  coupling  technique  of  the  two  methods  along  with  the 
difficulties  that  can  be  encountered  and  ways  of  circumventing  them. 

10.  Part  IV  gives  the  details  of  the  example  problem  solved  to  validate 
the  computer  code  developed.  The  results  are  presented  in  a  tabular  form  and 
a  few  cases  are  compared  in  the  complete  finite-element  version.  Discussion 
of  the  results  and  conclusions  is  presented  in  Part  V.  On  the  whole,  this 
procedure  is  found  to  be  very  effective  and  efficient  and  offers  greater  po¬ 
tential,  particularly  in  solving  hydraulic  U-lock  structure  problems.  The 
model  offers  capability  in  representing  large  linearly  elastic  soil  media  in 
solving  complex  soil-structure  interaction  problems  with  considerable  ease  in 
input-data  preparation.  Finally,  recommendations  for  future  work  are  listed. 


PART  II:  BOUNDARY-ELEMENT  METHOD 


Introduction 

11.  This  part  deals  with  the  basic  theory  of  the  BEM.  Equations  of 
linear  theory  of  elasticity  which  are  required  in  the  derivation  of  the 
boundary-element  equations  are  also  presented  briefly.  The  sequence  used  in 
the  derivation  is  for  the  direct  method  wherein  one  uses  a  fundamental  solu¬ 
tion  due  to  Kelvin  (Love  1927).  Other  simplified  boundary  element  techniques 
using  Flamant  equations  (Crouch  and  Starfield  1983)  are  also  investigated; 
details  of  this  method  were  gathered  in  a  preliminary  study  of  the  material 
herein.  The  boundary-element  equations  are  derived  in  such  a  form  that  they 
can  be  applied  to  plane  elasticity  problems  or  three-dimensional  elasticity 
problems.  The  difference  will  be  in  the  use  of  corresponding  fundamental 
solutions. 

Basic  Equations  of  Elasticity 

12.  The  basic  equations  of  elasticity  are  presented  here  for  a  three- 
dimensional  case.  Index  notation  is  used  for  convenience.  The  equilibrium 
equations  at  any  point  in  the  interior  domain  n  of  a  solid  continuum  are 


o.  .  .  +  b .  =  0  in  n 

ij ,  J  i 


where 

a . ,  =  stress  tensor 
b.  =  body  force  vector 

The  comma  after  ij  represents  differentiation  of  the  stress  tensor  o.^ 
with  respect  to  the  corresponding  axes  represented  by  the  subscript  following 
the  comma. 

13-  The  stress  tensor  has  to  match  with  the  prescribed  tractions  tL 
on  the  boundary,  i.e.. 


3 


where  n j  is  the  unit  normal  vector  on  the  boundary  .  Displacements  are 
prescribed  on  boundary  , 


u.  =  u.  on  r 
1  1 


(4 


In  this  case,  the  total  boundary  r  =  +  r2  .  Equations  2,  3,  and  4  become 

the  fundamental  equations  in  elasticity.  But  to  solve  these  equations,  one 
needs  two  additional  sets  of  equations:  one  set  to  represent  strain- 
displacement  relationships,  i.e., 


e .  .  =  x  (u.  .  +  u  .  . ) 

ij  2  i,j  J,i 


(5 


where  is  the  linear  strain  tensor,  and  the  second  set  to  represent  the 

stress  strain  relations,  i.e., 


O..  -  C  ,  (6 

ij  ijkl  kl 

where  is  a  fourth-order  material  property  tensor. 

Equations  for  Boundary-Element  Method 

14.  The  derivation  of  the  BEM  can  be  achieved  in  different  ways. 
Brebbia  and  Walker  (1972)  have  used  a  weighted  residual  approach  for  deriva¬ 
tion  of  the  boundary-element  equations.  The  approach  used  here  differs 

£ 

slightly  and  is  more  rational.  If  u.  is  assumed  as  a  weighting  function, 
then  the  equilibrium  equation  can  be  written  as 

f(o.  .  .  +  b. )  u*  dn  =  0  (7 

u  1J’J  1  1 

Considering  the  integral  shown  in  Equation  8, 

J(o.  ,u  )  .  dn  (8 

n  J 


and  differentiating  the  integrand  in  the  above  equation, 


lj  l.J 


J 


.  .  U  .  /  . 

lj  1  fJ 


uu 


J ' 
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iJ.J  1 


j  “ 
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By  using  the  divergence  theorem  and  substitution  in  Equation  3*  the  results 
are 


d!1 1  {(oijuI)nj 


dr 


=  Jt.u.  dr 
r  1  1 


The  substitution  of  Equation  10  into  Equation  9  and  rearranging  Equation  9 
into  the  form  of  equilibrium  Equation  7  results  in 


JVu*  dr  -  JojjU*  j  dn  ♦  J  t^u*  dn  =  o 


By  using  symmetry  of  stress  tensor,  i.e.,  o^  =  , 


dfl  !  K,j  *  “j.i’  d0 


=  Jo .  .e .  . 

n  1J 


dn 


The  substitution  of  Equation  12  into  Equation  7  results  in 


J  t.u*  dr  -  Jo.  e*  dn  +  Jb.u*  dn  =  0 
Jr  1  1  a  lj  lj  n  1  1 


By  using  the  stress-strain  relationship  of  Equation  6, 

hj'Ij d!1 "  £cuKie*ieij  d“ 


( 


Following  Equation  9, 


dQ  =  d!>  *  i0ijEij da  (,5) 


Following  Equation  10,  the  results  are 


J( 

n 


J 


dn 


Jo.  u.n 

r  ^  1  • 


dr 


(16) 


Substitution  of  Equation  16  into  Equation  15  and  using  Equation  14  gives 


Jo.  ,  ,u.  d£l  +  Jo. ,e .  , 
J  iJ.J  i  Ja  U  !j 


dn 


(17) 


Eliminating  the  second  integral  on  the  right-hand  side  of  Equation  17  by  using 
Equation  13  results  in 

Jt.u.  dr  -  Jt.u.  dr  +  Jo.,  .u.  dn  +  Jb.u.  dn  =  0  (18) 

r  r  1  n  n  1  1 

This  is  the  governing  equation  for  the  domain  under  consideration.  The  first 
domain  integral  in  the  above  equation  can  be  removed  by  assuming  a  solution  of 
an  equation  such  that 


°ij  j  +  Au  (s  ,  q)  =  0  in  o  (19) 


where 

s  =  source  point  where  a  unit  load  is  applied  in  the  1  direction 

q  =  field  point  where  the  displacements  and  tractions  are  calculated 
due  to  the  unit  load 

Mathematical ly , 

if  s  *  q  ,  Au  =  0  ; 

if  s  =  q  ,  and  i  =  1  ,  Ja ^  dn  =  1  , 

.  n  1 

and 

if  i  *  1  :  A . ,  =  0 


The  solutions  to  the  above  equation  u ^  and  t^  for  1  =  1  ,  2  ,  3 
represent  the  displacements  and  tractions  in  the  i  direction  due  to  a  unit 
concentrated  load  at  s  in  the  1  direction.  Substituting  Equation  19  into 
Equation  18,  for  any  source  point  s  ,  gives 

-uL(s)  +  Ju*ktk  dr  -  ftlkuk  dr  +  Julkbk  dn  =  0  (20 

r  r  a 

Omission  of  body  forces  leaves 

uL(s)  +  Jt*kuk  dr  =  Ju*ktk  dr  (21 

r  r 

Equation  21  is  for  a  source  point  inside  the  domain.  For  the  BEM,  the  source 

point  has  to  be  moved  on  to  the  boundary.  When  the  source  point  is  on  the 

boundary,  a  singularity  occurs  in  the  fundamental  solution  and  Equation  21 

has  to  be  integrated  in  a  special  manner.  In  Figure  3,  two  boundaries  are 

considered,  r^  for  r  =  e  at  the  source  point  and  r^  which  is  equal  to 

r  -  r  .  In  a  two-dimensional  case,  these  boundaries  are  lines,  and  for  a 
€ 


Figure  3.  Augmented  surface  for 
integration  on  the  boundary 

three-dimensional  case,  they  are  surfaces.  The  integration  has  to  be  per¬ 
formed  as  e  ■*  o  .  For  an  explanation  of  the  integration,  let  us  consider 


The  first  part  of  the  integral  in  Equation  22  can  be  shown  to  be  -  (u^/2)  , 
after  substitution  of  the  fundamental  traction  into  the  integral  sign  and 
noting  that  ear.  It  can  also  be  shown  that  the  second  part  does  not  in¬ 
troduce  any  new  term  when  £  -►  o  ,  when  Equation  21  is  used  on  the  boundary. 
Thus,  the  boundary  integral  equation  on  the  boundary  is  written  as 

csui(s)  +  Jt*kuk  dr  =  Ju*ktk  dr  (23) 

where  cg  =  1/2  when  boundary  is  smooth,  and  cg  =1  if  s  is  inside  the 
domain.  The  value  of  cg  for  higher  order  elements  can  be  derived  separately 
or  calculated  from  rigid  body  motion  criteria  (Brebbia  and  Walker  1972; 
Brebbia,  Telles,  and  Wrobel  1984). 

Fundamental  Solution 


15.  The  fundamental  solution  is  an  analytical  point-load  solution  in 
the  domain  and  thi~  is  used  to  convert  the  domain  integral  into  a  boundary 
integral.  The  solution  of  Equation  19  is  the  fundamental  solution  of  dis¬ 
placements  and  tractions  for  elastostatics  problems.  For  three-  and  two- 
dimensional  cases,  Kelvin  (Love  1927)  developed  the  fundamental  solution  due 
to  a  point  load  in  an  infinite  continuum.  Several  others  (Brebbia  and  Walker 
1972;  Brebbia,  Telles,  and  Wrobel  1984;  and  Hentenyi  1946)  produced  solutions 
for  different  domains  and  loads.  The  use  of  a  particular  fundamental  solu¬ 
tion  is  a  matter  of  choice  and  each  has  its  own  advantages  in  application. 
Kelvin's  solution  is  adopted  in  this  work  for  two  reasons:  first,  it  can  be 
used  for  bounded  domains;  second,  it  is  convenient  and  required  to  solve  prob¬ 
lems  in  layered  media.  The  expression  for  the  fundamental  solution  of  Kelvin 
for  displacements  and  tractions  for  a  unit  load  in  an  infinite  continuum  are 
given  below: 
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u.  .(s,q) 


_ 1 _ 

1 6 it  (1  -  v)  Gr 
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4  v )  6  .  .  +  r  .  r  > 

ij  .  i  »J 


(24) 


for  three-dimensional,  and 


^  .'TL'T 


uij<3’q)  =  s,<r.  v)c  {(3  - 4v>  in<r)4i  j  -  r,ir,j}  (2- 

for  two-dimensional  plane-strain  problems.  Corresponding  traction  functions 


t*  (s,q)  =  yj — - 7 —  if(1  -  2v)6 .  +  8r  .r  1 

lj  4an(1  -  v)ra  ij  .1  3n 


-  (1  -  2v)(r 


.n  -  r  ,n.  )> 

,i  j  ,J  i  f 


where  a  =  2  ,  1  ;  8=3,2  for  three-  and  two-dimensional  plane  strain, 

respectively.  Also,  r  =  r(s,q)  represents  the  distance  between  the  load 
point  s  and  the  field  point  q  and  its  derivatives  are  taken  with  reference 
to  the  coordinates  of  q  ,  i.e., 

r  =  (riri)1/2  =  | s  -  q |  , 

ri  =  xA(q)  -  x^s)  , 

„  r . 

3r  _  _i 

r,i  "  3x.(q)  "  r 

Interpolation  Functions  and  Numerical  Analysis 

16.  An  analytical  intergration  of  the  boundary-element  equation  as  seen 
in  Equation  23  is  extremely  tedious  and  not  practical  for  solving  engineering 
problems.  These  integrations  are  performed  in  a  piecewise  manner  on  the  boun¬ 
dary  using  discrete  boundary  elements.  The  boundary  is  discretized  into  line 
elements  in  a  two-dimensional  problem  and  surface  elements  in  a  three- 
dimensional  problem. 

17.  The  variations  of  the  unknown  displacements  and  tractions  on  the 
boundary  are  achieved  through  the  use  of  interpolation  functions,  namely, 

u.  =  ♦.  .U1? 
i  ij  J 

(27) 

t.  =  i>.  j" 

i  O  J 

where  U1^  and  T^  are  the  nodal  values  of  displacement  and  traction  vectors 


on  the  nv,‘  boundary  element,  respectively,  and  <j^  j  and  are  the 

interpolation  functions.  The  simplest  of  these  is  the  constant  element, 
where  the  displacements  and  tractions  are  considered  to  be  constants  within 
the  elements.  In  the  case  of  linear  elements,  the  displacements  and  trac¬ 
tions  vary  linearly  within  the  element.  Higher  order  elements  can  also  be 
formulated  (Banerjee  and  Butterfield  1981;  Brebbia,  Telles,  and  Wrcbel 
1984). 

18.  The  boundary  is  discretized  into  a  number  of  constant  elements. 
Then  Equation  23  is  applied  on  the  boundary  in  a  discrete  form.  The  corre¬ 
sponding  boundary-element  equation  in  the  constant  elements  would  be  of  the 
form, 

c(s)u  (s)  +  £  t*  dr\  U  =  £  u*  dr)  T  (2 

>'  b  )  (rj  ) 


where 


total  number  of  boundary  elements 
fch 

m  boundary  element 


rm  =  m  boundary  element 

Uj  and  Tj  =  displacements  and  tractions  in  the  element  j 

Equation  28  relates  the  ith  node  integration  (source  point  s)  and  the  in¬ 
tegration  is  carried  out  throughout  the  boundary  on  the  M  elements.  This 
is  done  sequentially  by  numerical  integration,  using  Gauss  quadrature  formu¬ 
las.  The  source  point  on  each  element  is  chosen  to  be  the  midpoint  of  the 
element.  Hence,  Equation  28  is  applied  at  every  element  at  its  midpoint  to 
form  the  system  equations.  The  integrals  over  each  element  now  become 
2x2  matrices  and  are  represented  as  H ■ j  and  G ^ j  .  Hence,  Equation  28 
becomes 


c. ( s ) U  (s) 


*  l  V,  =  I  GhTJ 

J=1  J  J  J=1  J  J 


Equation  29  relates  the  value  of  displacements  U  at  the  midpoint  of  the 
element  i  ,  namely  the  source  point  with  the  displacements  and  tractions  at 
all  the  elements  on  the  boundary,  including  the  source  point.  Equation  29 
can  be  written  as. 


where 


H. ,  =  H.  .  for  i  *  J 
lj  „ 1 j  J 

Hij  =  Hij  +  ci  for  1  =  J 

When  i  =  j  ,  in  Equation  29,  the  integration  becomes  singular  and  has  to  be 
evaluated  analytically  or  by  other  means.  It  is  easier  to  calculate  Ho¬ 
using  rigid  body  considerations  and  for  constant  element  it  works  out  to  be 
0.5  (Brebbia  and  Walker  1972).  can  be  calculated  analytically  or  by  a 

logarithmically  weighted  numerical  integration  formula. 

19.  The  system  equations  now  can  be  written  as 

[ H ] { U }  =  [ G  ]  { T }  (31) 

with  known  displacements  and  tractions,  specified  as  boundary  conditions. 

( H ]  and  [G]  are  square  matrices  of  the  order  2  x  M  . 

Solution  of  System  Equations 

20.  Whenever  displacements  are  prescribed,  the  tractions  cannot  be 
prescribed,  and  vice  versa.  Knowing  the  known  displacements  and  tractions, 
and  rearranging  the  set  of  equations  with  unknowns  on  the  left-hand  side  and 
knowns  at  the  right-hand  side,  the  entire  Equation  31  is  rearranged  as  a  set 
of  linear  simultaneous  equations. 


( A ]  { x )  =  fb]  (32) 

Solving  Equation  32,  we  get  all  the  unknown  displacements  and  tractions  on  the 
boundary . 


PART  III:  COUPLING  OF  BOUNDARY-  AND 
FINITE-ELEMENT  METHODS 


Introduction 


21.  The  BEM  has  been  found  to  be  a  very  powerful  technique  for  solv¬ 
ing  stresses  and  deformations  of  large  continuum  where  the  boundaries  are  at 
large  distances.  This  method  is  used  to  represent  large  soil  media  in  the 
soil-structure  interaction  problems.  The  FEM  is  used  to  represent  the  struc¬ 
ture;  properties  such  as  complex  geometry,  heterogeneity,  nonlinearity,  rein¬ 
forcement,  etc.,  are  better  modeled  by  the  FEM.  For  an  elastic  continuum, 
the  BEM  yields  better  results  with  relatively  lesser  number  of  unknowns  with 
easy  preparation  of  input  data.  If  one  needs  to  use  nonlinear  properties  of 
the  soil  in  the  immediate  vicinity  of  the  structure,  it  is  recommended  at 
present  to  use  the  FEM  for  the  nonlinear  portion  of  the  soil  media.  The  use 
of  these  two  methods  make  the  analysis  computationally  and  economically  more 
efficient . 


Equations  Used  in  Coupling 

22.  There  are  two  basic  procedures  for  combining  the  two  methods 
(Brebbia  and  Walker  1972;  Brebbia,  Telles,  and  Wrobel  1984).  One  is  to 
convert  the  finite-element  equations  to  suit  the  boundary-element  equations. 
This  has  been  the  main  technique  used  by  many  researchers.  However,  this 
procedure  has  serious  efficiency  problems,  especially  when  the  FEM  has  non¬ 
linear  behavior.  The  method  adopted  here  is  in  the  reverse  order;  i.e., 
the  boundary-element  equation  is  transformed  into  an  equivalent-stiffness 
matrix  and  added  to  the  finite-element  half-banded  stiffness-matr i x  equa¬ 
tion.  This  technique  has  many  advantages  in  solving  so i 1 -structure  inter¬ 
action  problems. 

23.  Consider  two  regions,  and  ,  with  different  material  prop¬ 
erties  and  geometries  as  shown  in  Figure  4.  Region  is  divided  into 

finite  elements  and  region  into  boundary  elements.  Coupling  the  two 

models  is  done  by  using  equilibrium  and  compatibility  on  the  interface  : 

a.  Equilibrium  of  tractions;  i.e.,  the  tractions  on  inter¬ 

face  for  region  ft  and  q  should  be  equal  and  opposite. 


v  BEM  Model 

Figure  4.  Domain  with  different  properties 


b.  Compatibility  of  displacements;  i.e.,  displacements  on  gg 
interface  for  region  and  should  be  equal. 

Finite-Element  Equation 

24.  The  structure  is  discretized  into  plane-strain  finite  elements  in 
analysis.  The  structure  stiffness  matrix  is  developed  using  the  compute 
,  and  this  includes  the  soil-structure  interface,  denoted  by  gg  in  Fig 
4.  The  finite-element  stiffness-matrix  equation  will  be 


soil  stiffness 


KSS  = 

KSB’  KBS  =  "diagonal  stiffness  beams 
Kgg  =  beam  stiffness 
Ug  =  soil  displacements 
Ug  =  beam  displacements 
Fg  =  forces  on  soil 
Fg  =  forces  in  beam 

where  the  square  matrix  is  the  global-stiffness  matrix  of  the  finite-element 
portion,  Ug  and  Fg  are  the  displacement  and  force  vectors  of  the  boundary- 
element  interface  and  Ug  and  Fg  are  the  remaining  displacement  and  force 
vectors  of  the  finite-element  system.  The  finite-element  stiifness  matrix 
shown  in  Equation  33  is  developed  as  a  half-banded  matrix  for  computational 
economy . 


Boundary-Element  Equation 

25.  Using  the  BEM,  a  stiffness  matrix  representing  the  soil  media  has 
to  be  developed  on  the  interface.  This  should  be  in  such  a  form  that  it  is 
compatible  with  the  finite-element  system.  Therefore,  the  boundary-element 
stiffness-matrix  equation  should  be  of  the  type: 

[kg] {Ug}  =  -{Fg}  ♦  {fB}  (34) 

where  fg  is  a  condensed  force  vector  representing  the  prescribed  forces  and 
displacements  on  the  boundary  of  the  soil  medium.  The  development  of  the 
condensed  equation  (Equation  34)  representing  the  soil  medium  for  a  hydraulic 
U-lock  structure  is  discussed  in  the  following  section. 

Soil  Stiffness 


26.  To  obtain  a  condensed  stiffness  matrix  [kg]  to  represent  soil, 
Kelvin's  fundamental  solution  as  shown  in  Equations  24  through  26  is  used  for 
the  formulation  of  the  boundary-element  equations.  Figure  5  shows  the  BEM  for 
a  U-lock  structure.  Advantage  in  symmetry  of  the  problem  is  taken  care  of  in 
the  computer  code. 


Figure  5.  Boundary-element  model  for  a  U-lock  structure 

27.  The  boundary-element  equations  are  developed  using  Equation  26  and 
the  Final  set  of  matrices  will  be  of  the  form: 


H-Matrix  G-Matrix 


where  U  and  T  are  the  nodal  d  isplacemer.t  and  traction  vectors  on  the  four 
boundaries  gg  ,  g1  ,  g2  ,  and  as  shown  in  Figure  5.  The  subscripts 

on  U  and  T  denote  the  respective  boundaries.  The  following  boundary  con¬ 
ditions  are  prescribed  for  the  problem: 
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’  T1x 
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,  a  prescribed  traction; 
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'  U2x 

=  0 

• 

r2v 

-  0 

* 

on 
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'  u3v 
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T3x 
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;  and 

on 

fa  . 

,  both 

Ur 

and  T„  are  unknowns. 

Considering  the  prescribed  and  the  unknown  components  of  displacements  and 
tractions,  the  H  and  G  matrices  can  be  reorganized  such  that  the  new  H 
and  G  matrices  are: 


(35) 


where  Vp  and  Vy  are  vectors  representing  the  prescribed  and  the  unknown 
components  of  all  displacements  and  tractions  on  g i  f  g2  ,  and  g^  .  The 
sizes  of  submatrices  of  H  and  G  correspond  to  the  orders  of  the  interface 
vectors  Ug  and  Vy  ,  respectively.  From  Equation  35 


Vy  =  H22(G21TB  +  G22VP  *  H21UB 


(36) 


and  substituting  this  into  part  of  Equation  35, 


(H 


11 


H12H22H21)UB 


(G 


1 1 


-  h12H22G21)TB  +  (G 


12 


-  H12H22G22)VP  (37; 


this  equation  can  be  reduced  to 


[kg]{UB)  =  -{Fg}  +  {fg}  (38) 

where  kg  is  the  required  stiffness  matrix  of  the  soil  medium,  Fg  is  the 
equivalent  nodal  forces  of  the  traction  vector  Tg  ,  and  fg  is  an  equivalent 
force  vector  representing  the  prescribed  forces  and  displacements  on  gi  , 
g2  ,  and  g^  .  The  development  kg  and  fg  appears  to  be  quite  cumbersome. 
But,  the  computations  can  be  simplified  considerably  by  using  a  static  conden¬ 
sation  procedure  (Vesic  1961).  Using  the  Gauss/Jordan  elimination  procedure, 
if  one  transforms  both  H  and  G  matrices  such  that  H22  is  made  into  an 
identity  matrix  and  H^2  to  a  null  matrix,  then  the  new  transformed  matrices 
H.|j  ,  G^  ,  and  G^2  are,  i.e., 

H1 1UB  =  G11TB  +  G12VP  ( 

the  matrices  within  the  respective  brackets  in  Equation  36.  Again  using 


23 


is  made  into 


a 


elimination  procedure,  Equation  39  is  transformed  such  that 
an  identity  matrix,  i.e., 


h11uB  =  tB  +  G12VP 


Transforming  nodal  tractions  vector  Tg  to  equivalent  nodal  force  vector  -Fg 
in  Equation  40  and  adding  this  to  Equation  33,  we  get  the  desired  Equation  42. 
This  is  achieved  by  multiplying  both  sides  of  Equation  40  by  a  distribution 
matrix  M  ,  to  convert  tractions  into  nodal  forces.  Thus,  we  have 


[MHn]UB  =  MTg  +  MG12Vp 


-  ~{fn}  +  {fR} 


i.e. , 


[kB]{UB}  = 


(34  bis) 


Coupling  of  Finite-Element  and  Boundary-Element  Matrices 

28.  Having  developed  the  stiffness  matrices  of  the  structure  and  the 
soil  as  shown  in  Equations  33  and  34,  the  coupling  is  done  ensuring  compati¬ 
bility  and  equilibrium  at  the  interface.  Now  the  combined  FEM/BEM  model  can 
be  written  as 


ss 

CO 

c n 

BS 

kbb 

This  equation  is  readily  solved  to  get  the  displacements  and  stresses  of  the 
structure. 

Asymmetry  of  Boundary-Element  Stiffness  Matrix 

29.  The  stiffness  matrix  [kg]  formed  from  the  boundary-element  domain 
is  generally  asymmetric.  This  poses  problems  for  addition  into  the  symmetric¬ 
stiffness  matrix  of  the  finite-element  region.  The  asymmetry  of  the  boundary 
element  is  not  new  and  has  been  reported  in  literature  before  (Banerjee  and 
Butterfield  1981;  Brebbia,  Telles,  and  Wrobel  1984;  Georgiou  1981;  Hartmann 
1981).  This  has  been  attributed  to  three  factors;  namely,  discretization  of 
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s  s 


the  boundary-element  domain,  collocation  process,  and  the  nature  of  the  fun¬ 
damental  solution.  To  alleviate  this  anomaly,  the  boundary-element  stiffness 
matrix  is  symmetrized,  discarding  the  unsymmetric  part.  This  has  been  justi¬ 
fied  by  Georgiou  (1981)  by  showing  examples,  the  results  of  which  are  reason¬ 
ably  accurate.  A  detailed  mathematical  procedure  to  derive  a  symmetric  stiff¬ 
ness  matrix  is  explained  by  Hartmann  (1981)  in  that  by  solving  G-1Hu  =  t 
with  Galerkin's  method  or  by  minimizing  a  potential  <t>(u)  ,  a  symmetric 
stiffness  matrix  is  guaranteed. 

30.  In  this  present  work,  the  boundary-element  stiffness  matrix  is 
symmetrized  by  averaging  the  off-diagonal  terms,  by  the  principle  of  least 
squares  as  suggested  by  Brebbia  and  Walker  (1972).  It  is  found  that  the  error 
involved  is  very  small  and  negligible  that 

k  .  .  -  k  . . 

^  <  0.002£ 
i  i 


Compatible  Boundary  Stiffness  Matrix 

31.  The  boundary  elements  representing  the  soil  are  assumed  to  have  a 
constant  variation  of  displacements  and  traction  over  their  length.  Thus,  on 
the  interface  the  boundary  elements  will  lie  in  between  the  finite-element 
nodes.  The  condensed  effective  stiffness  of  the  soil  portion  obtained  by 
Equation  39  has  to  be  made  compatible  with  the  finite-element  stiffness  matrix 
representing  the  structure  for  coupling. 

Coupling  of  stiffness  affected 

32.  The  coupling  of  stiffness  from  the  constant  boundary-element  form¬ 
ulation  of  the  soil  to  the  finite-element  stiffness  of  the  structure  is  af¬ 
fected  in  two  ways,  as  explained  below. 

a.  Method  I .  In  this  method,  the  midpoints  of  the  boundary  ele¬ 
ments  are  positioned  in  such  a  way  that  they  lie  on  the  finite- 
element  nodes  of  the  structure  at  the  interface.  The  stiffness 
at  the  ends  of  the  boundary  elements  at  the  interface  are  re¬ 
duced  by  half  so  as  to  correspond  to  the  stiffness  of  the  sub¬ 
grade  at  the  interface,  to  be  added  appropriately  at  the 
interface . 

b.  Method  1 1 .  In  this  method,  the  end  points  of  the  boundary 
elements  coincide  with  the  finite-element  nodes.  Therefore, 
the  stiffness  matrix  of  the  subgrade  from  the  boundary-element 
region  will  be  lesser  by  two  rows  and  two  columns  due  to  a  node 


Subgrade 
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lesser  in  boundary-element  formulation.  The  condensed  stiff¬ 
ness  matrix  at  the  interface,  representing  the  subgrade,  is 
expanded  by  premultiplying  it  by  a  transformation  matrix  and 
then  postmultiplying  it  by  its  transpose.  This  transformation 
is  possible  by  the  use  of  contragradient  law  (Vallabhan  and 

Jain  1972)  and  satisfies  energy  principles.  Let  Ub  and  Fb 
be  the  nodal  displacements  and  forces  at  the  midpoints  of  the 

f  f 

boundary  elements  at  the  interface.  Let  U  and  F  be  the 
displacements  and  forces  at  the  finite-element  nodes  at  the  in¬ 
terface.  Let  [T]  be  the  transformation  matrix  which  converts 
the  finite-element  displacements  of  the  interface  nodes  into 
the  boundary-element  displacement  at  the  midpoints.  Therefore, 

{ub}  =  [T]{uf}  (43) 


By  contragradient  law, 

{ Ff }  =  [ T ] 1  ( Fb}  (44) 


Transposing  Equation  46  and  postmultiplying  by  { u ^ )  , 

{FfMuf}  =  {FbMT]{uf}  =  {Fb}fc{ub)  (45) 


i.e.,  the  work  done  by  the  forces  on  the  displacements  of  the 
finite-element  systems  is  equal  to  the  work  done  by  the  forces 
on  the  displacements  of  the  boundary-element  system. 


matrix 

The  condensed  stiffness  matrix  of  the  subgrade  is 


[kB]{UB}  -  fB 
[ kB ] [T ] [UF]  =  {Fb} 

[ TT ] [ kB ] [T ] (UF }  ,  [TT]{Fb}  =  { Ff} 
[ kR ] { Uf }  =  { FF } 


(46) 


This  method  is  more  convenient  and  yields  a  better  solution  than  Method  I. 


PART  IV:  SOIL-STRUCTURE  INTERACTION  PROBLEMS 


A  Typical  Problem 

34.  Many  researchers  have  validated  the  use  of  the  BEM  for  solving 
stresses  and  displacements  in  linear  continuum  (Banerjee  and  Butterfield  1981; 
Brebbia  and  Walker  1972;  and  Brebbia,  Telles,  and  Wrobel  1984)  and  hence,  ex¬ 
amples  are  not  shown  here  for  demonstrating  the  capabilities  of  the  boundary- 
element  method.  On  the  other  hand,  examples  are  shown  here  to  illustrate  the 
coupling  of  the  BEM  and  FEM  for  solving  soil-structure  interaction  problems. 

A  typical  problem  which  concerns  a  hydraulic  U-lock  structure  was  selected  for 
this  purpose.  This  particular  problem  has  been  studied  and  solved  by  fellow 
researcher,  Wilson  (1974),  who  used  other  methods  of  analyses.  The  problem 

deals  with  a  strip  of  plate  of  unit  width  resting  on  finite-soil  media  and  is 

considered  as  one  of  place-strain  type.  The  strip  is  modeled  as  a  beam  of 
length,  L  equal  to  100  ft.*  Two  cases  of  depth  d  of  beam  are  considered, 
one  for  depth  equal  to  5  ft  and  the  other  equal  to  10  ft.  Two  values  of  the 
depth  of  the  soil,  H  ,  were  used;  in  the  first  case,  H  is  set  equal  to  100 
ft,  and  in  the  second  case,  H  is  set  as  200  ft.  The  modulus  of  elasticity 
of  the  concrete  beam  Eb  is  taken  as  3  *  10^  psi  .  The  modulus  of  elastic¬ 
ity  of  the  soil  Es  is  varied  such  that  the  two  values  of  the  modular  ratio 

n  =  Eb/Eg  are  taken  as  10  and  100,  respectively.  For  convenience,  in  this 

analysis  the  Poisson's  ratio  of  the  concrete  and  the  soil  is  set  equal  to  0.2, 
even  though  this  is  not  a  limitation  of  the  methodology  at  all.  Altogether, 
there  are  eight  cases  of  beam/soil  geometry  and  material  properties  for  each 
loading  case.  Table  1  indicates  the  nomenclature  for  these  eight  cases.  The 
load  cases  used  here  are: 

a.  A  uniformly  distributed  load  of  5  k  ft  on  top  of  the  he  m . 

b.  A  concentrated  load  of  250  k  at  each  j  oj'  tr.e  c*  «r  . 

c.  A  concentrated  moment  of  5,0('0  ft-k  i*  ...  <  ;  • 

It  is  further  assumed  that  the  soil  io  resting  or  1  snoot  r  , 

surface.  The  effect  of  the  smoothness  of  the  hot:.."’  :>,•  lev  •  *  he  .ecu:  . -v 
of  the  results  is  found  to  be  negligible  for  the  ,r.  ...  c-m-e;,.  .  r  c>»  w'  ire 


*  A  table  of  factors  for  converting  non -SI  ur.  1  *  s  of  mo  i.;:,r‘'t>'t.t  to  SI 
(metric)  units  is  presented  on  pngo  1. 


dealing  with  finite  dimensions  of  soil  media,  it  is  assumed  that  the  vertical 
boundaries  at  a  distance  of  200  ft  away  from  the  center  line  are  considered 
smooth.  All  these  boundary  conditions  are  selected  only  from  a  practical  de¬ 
sign  point  of  view,  and  they  reflect  no  limitations  of  the  technique.  The 
loads,  geometry,  and  the  boundary  conditions  of  the  problem  are  illustrated 
in  Figure  6.  Thus  far,  we  have  24  problems  to  be  solved. 

250  k  250  k 


E5  =  3,000  ksi 

vjj  =  vs  =  0.2 

H  =  100  ft,  200  ft 


L  =  100  ft 
L/.  =  10,  20 


-/t 


Eb 


n  =  r-  «  10,  100 

Bottom  and  sides  are  smooth  boundaries 


Figure  6.  Details  of  the  example  problem 


35.  Another  question  which  comes  into  the  analysis  is  the  compatibility 
of  horizontal  displacements  between  the  structure  and  the  soil  at  the  inter¬ 
face.  It  is  a  common  practice  in  beam-on-elastic-foundation  studies  to  ignore 
the  compatibility  of  horizontal  displacements  at  the  interface.  If  one  ig¬ 
nores  this  condition,  we  are  essentially  assuming  that  the  soil  interface  is 
perfectly  smooth.  Thus,  the  analysis  is  performed  in  two  categories.  In  the 
first  category,  the  compatibility  of  horizontal  displacements  are  neglected, 
thereby  making  the  soil  as  a  smooth  boundary.  In  the  second  case,  compati¬ 
bility  of  horizontal  displacements  on  the  interface  are  enforced. 


Discretization  of  the  Example  Problem 


36.  Due  to  the  symmetry  of  the  problem  only  half  the  domain  was  used 
for  discretizations.  Rectangular  finite  elements  with  8  deg  of  freedom  are 
employed  here.  The  beam  is  discretized  into  five  layers  with  ten  elements  in 
each  layer.  This  discretization  is  found  to  model  the  accurate  bending  of 
beams  for  displacements  and  stresses.  The  boundary-element  discretization 
must  be  made  compatible  with  the  finite-element  discretization  at  the  inter¬ 
face.  The  sizes  of  the  boundary  elements  are  varied  such  that  wherever  the 
displacements  and  stresses  are  small  or  uniform,  larger  elements  are  used  as 
shown  in  Figure  7.  A  convergence  study  was  made  on  the  number  of  boundary 
elements  in  the  combined  model,  and  it  was  found  that  the  discretization  with 
40  elements  shown  in  Figure  7  gave  essentially  the  same  results  with  larger 
numbers  of  boundary  elements. 

37.  In  order  to  verify  the  accuracy  of  the  finite-element/boundary- 
element  model,  a  complete  finite-element  study  is  also  made  for  four  selected 
cases,  such  as  5,  6,  7,  and  8  for  each  loading  condition.  The  finite-element 
discretization  for  the  problem  in  the  loadings  and  boundary  conditions  are 
shown  in  Figure  8. 


Presentation  of  Results 


38.  The  results  for  various  cases  are  presented  in  graphical  and  tabu¬ 
lar  form.  They  relate  principally  to  the  vertical  displacements  and  tractions 
at  the  interface. 

39-  These  results  are  seen  in  a  clearer  perspective  in  graphical  form. 
Parts  a  of  Figures  9  through  20  show  the  plots  of  vertical  displacement  at  the 
interface  for  each  loading  condition.  Parts  b  of  Figures  9  through  20  show 
the  plots  of  corresponding  vertical  tractions  at  the  interface.  These  graphs 
are  analyzed  for  the  first  category,  i.e.,  where  the  compatibility  of  horizon¬ 
tal  displacements  are  not  enforced  at  the  interface. 

40.  In  all  tables,  the  displacements  shown  at  nodes  1  to  1 1  are  from 
the  line  of  symmetry  at  5-ft  intervals,  at  the  soil-structure  interface. 

In  the  case  of  vertical  traction  at  the  soil-structure  face,  the  values  are 
a  constant  in  each  of  the  boundary  elements  and  are  shown  at  the  midpoint 
of  each  boundary  element.  Tables  A - 1 ,  A-3,  and  A-5  display  the  vertical 
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Figure  9a.  Vertical  displacements  at  the  interface  for  cases  1  and  2 
with  uniformly  distributed  load  5  k/ft  on  top  of  beam 


Figure  9b.  Vertical  tractions  at  the  interface  for  cases  1  and  2 
with  uniformly  distributed  load  5  k/ft  on  top  of  beam 
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Figure  12a.  Vertical  displacements  at  the  interface  for  cases  7  and  8 
with  uniformly  distributed  load  5  k/ft  on  top  of  beam 


O'- CASE  I 
•  -CASE  2 


Case  1 

1 

10 

10 

Case  2 

1 

10 

100 

FE/BE  Interface 

- >  x  (ft) 


Figure  13a.  Vertical  displacements  at  the  interface  for  cases  1  and  2 
with  concentrated  load  of  250  k  at  each  end  of  the  beam 
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Figure  13b.  Vertical  tractions  at  the  interface  for  cases  1  and  2 
with  concentrated  load  of  250  k  at  each  end  of  the  beam 


Vertical  Tractions,  k/fl 


Figure  14b.  Vertical  tractions  at  the  interface  for  cases  3  and  4 
with  concentrated  load  of  250  k  at  each  end  of  the  beam 
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Figure  15a.  Vertical  displacements  at  the  interface  for  cases  5  and  6 
with  concentrated  load  of  250  k  at  each  end  of  the  beam 
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Figure  15b.  Vertical  tractions  at  the  interface  for  cases  5  and  6 
with  concentrated  load  of  250  k  at  each  end  of  the  beam 
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Figure  16b.  Vertical  tractions  at  the  interface  for  cases  7  and  8 
with  concentrated  load  of  250  k  at  each  end  of  the  beam 
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Figure  17a.  Vertical  displacements  at  the  interface  for  cases  1  and  2 
with  concentrated  moment  of  5,000  ft-k  at  end  end  of  beam 
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displacements  at  the  interface  and  Tables  A-2,  A-4,  and  A-6  display  the  ver¬ 
tical  tractions  at  the  interface  for  the  first  category.  It  is  in  this  cate¬ 
gory  that  the  horizontal  displacements  at  the  soil  interface  are  neglected. 
Similarly,  vertical  displacements  and  tractions  at  the  interface  for  the  sec¬ 
ond  category  are  listed  in  Tables  B-1  through  B-6.  The  compatibility  of  the 
vertical  and  horizontal  displacements  of  the  soil  at  the  interface  are  en¬ 
forced  in  this  category.  The  results  of  the  complete  finite-element  study  for 
cases  5,  6,  7,  and  8  for  H/L  =  2  for  each  loading  condition  are  shown  in 
Tables  C-1  through  C-3.  Here,  only  the  displacements  of  the  soil-structure 
interface  are  shown. 

41.  The  comparisons  of  interface  displacements  obtained  from  the 
finite-element  analysis  for  cases  5,  6,  7,  and  8  for  H/L  =  2  for  each 
loading  condition  are  made  with  the  corresponding  cases  in  the  coupled 
FE/BE  model  of  the  second  category.  The  horizontal  displacements  of  the  soil 
interface  are  also  included  in  the  analysis  and  the  plots  are  shown  in 
Figures  21  through  26. 
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PART  V:  SUMMARY  AND  CONCLUSIONS 


Summary  and  Findings 

42.  The  feasibility  study  of  the  boundary-element  technique  for  stress 
analysis  of  soil  media  was  made  by  Vallabhan  (1983).  Based  on  promising  re¬ 
sults  for  elastostatic  problems,  the  research  was  continued  for  application 
specifically  in  soil-structure  interaction  problems.  Here,  the  structure  is 
modeled  by  the  BEM  and  the  soil  is  modeled  by  the  FEM.  Since  the  soil  is  con¬ 
sidered  as  a  two-dimensional  plane-strain  problem,  the  boundary  consists  of 
lines  and  input  data  for  these  problems  are  extremely  simple.  A  model  problem 
suggested  by  the  Waterways  Experiment  Station  has  been  solved.  Results  are 
compared  with  full  finite-element  solutions.  The  comparisons  are  very  good. 

43.  Since  the  methodology  and  the  equations  obtained  from  these  two 
numerical  methods  are  different,  the  emphasis  on  this  research  was  to  develop 
a  convenient  and  economical  technique  to  couple  the  two  methods  of  analysis 
for  solving  soil-structure  interaction  problems.  Using  a  static  condensation 
procedure,  the  boundary-element  equations  are  reduced  and  transformed  into  a 
stiffness  matrix  on  the  interface  boundary  between  the  structure  and  the  soil 
medium.  This  way  the  properties  of  the  half-band  width  of  the  finite-element 
stiffness  matrix  is  preserved.  A  computer  code  for  the  coupling  of  the  two 
methods  has  been  developed  in  FORTRAN. 

Conclusions 

44.  The  following  conclusions  are  made  from  this  study: 

a.  The  results  of  the  FEM/BE  combined  model  compare  very  well  with 
the  full  FEM  model  for  the  soil  and  the  structure. 

b.  The  coupling  method  developed  here  is  found  to  be  very  efficent 
in  numerical  computation  and  in  preparation  of  input  data. 

c.  Constant  boundary  elements  are  employed  here.  The  stiffness 
matrix  computed  using  constant  boundary  elements  is  almost 
symmetric  and  the  maximum  error  in  lack  of  symmetry  when  com¬ 
pared  to  the  diagonal  elements  for  the  above  problem  is  less 
than  0.002  percent. 
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Table  1 

Nomenclature  for  Eight  Loading  Cases  of  Beam/Soil 
Geometry  and  Material  Properties 


\\S\ 


Table  A1 

Vertical  DisDlacements  at  Interface,  Loading  Case  1 


Table  A3 

Vertical  Displacements  at  Interface,  Loadinc  Case  2 
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=  10 

L/t 

=  20 

L/t 

=  10 

L/t 

=  20 

Cx  £1  1  X 

n  =  i 

Eb/Es 

n  r  , 

Eb/Es 

n  =  1 

VEs 

n  -  , 

ek/e_ 

D  S 

10 

100 

10 

100 

10 

100 

10 

100 

Case 

1 

2 

3 

4 

5 

6 

7 

8 

xIO'3 

,io-2 

x10‘3 

x10-2 

x10‘3 

,10‘2 

«10‘3 

,  io-2 

Node  No. 

1 

-2.924 

-4.748 

-3.328 

-2.953 

-6.616 

-8.443 

-7.028 

-6.644 

2 

-2.983 

-4.817 

-3.354 

-3.028 

-6.674 

-8.511 

-7.051 

-6.717 

3 

-3.169 

-5.024 

-3.430 

-3.258 

-6.854 

-8.714 

-7. 120 

-6.941 

4 

-3.508 

-5.368 

-3.564 

-3.663 

-7.182 

-9.053 

-7.242 

-7.337 

5 

-4.043 

-5.850 

-3.776 

-4.275 

-7.703 

-9.527 

-7.439 

-7.936 

6 

-4.835 

-6.466 

-4.130 

-5.135 

-8.478 

-10.130 

-7.773 

-8.780 

7 

-5.964 

-7.214 

-4.770 

-6.29 2 

-9.588 

-10.870 

-8.392 

-9.921 

8 

-7.523 

-8.087 

-5.984 

-7.795 

-11.130 

-11.730 

-9.583 

-11.410 

9 

-9.636 

-9.076 

-8.293 

-9-681 

-'3.220 

-12.710 

- ’ ' . 8”0 

-  13.270 

10 

-12.230 

-10.150 

-12.310 

-11.950 

-15.790 

-13.780 

-15.870 

-15.520 

1 1 

-15.090 

-11.260 

-18.320 

-14.490 

-18.640 

-14.870 

-21.870 

-18.070 

Table  A5 

Vertical  Displacements  at  Interface,  Loading  Case  3 


Moment  of  5,000  ft-k — Beam  Ends 


Note:  L  ;  100  ft  ,  Eh  =  432,000  k/ft2  ,  v.  =  0.2  ,  v  ;  0.2 


APPENDIX  B:  VERTICAL  DISPLACEMENTS  AND  TRACTIONS 
AT  INTERFACE,  CATEGORY  2 
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Table  B1 

Vertical  Displacements  at  Interface.  Loading  Case  1 
UDL  5  k/ft — Top  of  Beam 


Note:  L  =  100  ft  ,  Eh  =  432,000  k/ft2  ,  v.  =  0.2  ,  v  =  0.2 


& 


Table  B3 

Vertical  Displacements  at  Interface,  Loading  Case  2 


Point  Load  250  k--Beam  Ends 


Note:  L  =  100  ft  ,  Eb  ;  U32.000  k/ft2  ,  =  0.2  ,  vg  =  0.2  . 


BA 


Element  No. 


1 

0.746 

-0.197 

0.062 

1 . 150 

0.818 

-0.071 

0.078 

2 

0.768 

-0.328 

0.109 

1.130 

0.843 

-0.206 

0.126 

3 

0.792 

-0.600 

0.227 

1 .050 

0.871 

-0.488 

0.247 

4 

0.764 

-1.040 

0.465 

0.847 

-0.847 

-0.940 

0.491 

5 

0.578 

-1.680 

0.886 

0.408 

0.664 

-1.610 

0.921 

6 

0.042 

-2.600 

1.490 

-0.464 

0. 126 

-2.560 

1.530 

7 

-1.200 

-3-920 

2.030 

-2.100 

-1.140 

-3.920 

2.  100 

8 

-4.260 

-5.980 

1 .090 

-5.160 

-4.250 

-6.050 

1 .  160 

9 

-10. 100 

-9.090 

-4.660 

-10.400 

-10.200 

-9.280 

-4.660 

10 

-38. 100 

-24.600 

-51.700 

-36.400 

-38.600 

-24.900 

-52.000 

Table  B5 

Vertical  Displacements  at  Interface.  Loading  Case  3 
Moment  of  5,000  ft-k--Beam  Ends 


Node  No. 


1 

-3-371 

-2.809 

-2.684 

-3-990 

-3.246 

-2.769 

-2.567 

2 

-3.391 

-2.758 

-2.733 

-3.997 

-3.266 

-2.719 

-2.616 

3 

-3.440 

-2.605 

-2.893 

-4.004 

-3.316 

-2.565 

-2.777 

4 

-3.487 

-2.340 

-3.200 

-3-970 

-3.364 

-2.301 

-3.084 

5 

-3.470 

-1.953 

-3.704 

-3.820 

-3-349 

-1.914 

-3.590 

6 

-3.284 

-1.423 

-4.438 

-3.440 

-3.164 

-1.385 

-4.325 

7 

-2.766 

-0.732 

-5-288 

-2.661 

-2.648 

-0.694 

-5. 177 

8 

-1 .657 

0.148 

-5.751 

-1.257 

-1.541 

0. 186 

-5.643 

9 

0.658 

1  .270 

-4.267 

1 .066 

0.773 

1.307 

-4.161 

10 

4.950 

2.691 

3.580 

4.751 

5.064 

2.729 

3.685 

1 1 

10.660 

4.292 

22.890 

9.937 

10.780 

4.329 

22.990 

Note:  L  =  100  ft  ,  Eu  =  432,000  k/ft^  ,  --  0.2  ,  v  =  0.2 


Table  Cl 

Vertical  Displacements  at  Interface,  Loading  Case  1 
UDL  5  k/ft — Top  of  Beam 


Table  C2 

Vertical  Displacements  at  Interface,  Loading  Case  3 
Moment  of  5,000  ft-k — Beam  Ends 


Body  force  vector 

^ijkl  Fourth-order  material  property  tensor 

432,000  k/ft^,  modulus  of  elasticity  of  beam 

fg  Condensed  force  vector  representing  prescribed  forces  and 
displacements  on  the  boundary  of  the  soil  medium 

Fg  Forces  in  beam 

Fs  Forces  on  soil 

[G]  and  [H]  Square  matrices  of  the  order  2  x  m 

k  Proportionality  constant,  known  as  the  subgrade  modulus  or  the 

modulus  of  subgrade  reaction 

[kB]  Condensed  stiffness  matrix  to  represent  soil 

KgB  Beam  stiffness 

Kgg,  Kgg  Off~diagonal  stiffness  beams 

Kgg  Soil  stiffness 

L  Length  of  beam 

M  Total  number  of  boundary  elements 

n  Unit  normal  vector  on  the  boundary 

p  Pressure  acting  on  the  soil  surface 

q  Field  point  where  the  displacements  and  tractions  are 
calculated  due  to  the  unit  load 

r  =  r(s,q)  Distance  between  the  load  point  s  and  field  point  q 

s  Source  point  where  a  unit  load  is  applied  in  the  specified 

direction 

t^  Prescribed  tractions  on  boundary 

» 

tii  Tractions  in  the  i  direction  due  to  a  unit  concentrated  load 
at  s  in  the  1  direction 

[T]  Transformation  matrix  which  converts  finite-element 

displacements  of  interface  nodes  into  the  boundary-element 
displacement  at  the  midpoints 

T^  Tractions  in  element  i 

i.  L. 

Tj  Nodal  values  of  traction  vectors  on  the  nun  boundary  element 

Tg  Traction  vector 

U|  Assumed  weighting  function 

* 

u^  Displacements  in  the  l  direction  due  to  a  unit  concentrated 
load  at  s  in  the  1  direction 

Ub  and  Fb  Nodal  displacements  at  midpoints  of  the  boundary  elements  at 
the  interface 


VD  and  V, 


,  and  i>.  , 
J 


Displacement  in  element  j 

Nodal  values  of  displacement  on  the  nth  boundary  element 
Beam  displacements 
Soil  displacements 

Vectors  representing  the  prescribed  and  the  unknown  components 
of  all  displacements  and  tractions  on  boundaries  g1  ,  g2  , 
and  g3 

Deflection  of  the  loaded  region  on  the  surface 
Complete  exterior  boundary 
Boundary  for  element 
Exterior  boundary  except  for 

Circular  augmented  surface  at  source  point  for  radius  (r) 
equal  to  some  assumed  error  (e) 

Exterior  boundary  where  tractions  are  prescribed 

Exterior  boundary  where  displacements  are  prescribed 

Linear  strain  tensor 

Stress  tensor 

0.2,  Poisson's  ratio  of  beam 
0.2,  Poisson's  ratio  of  soil 
Interpolation  functions 
Interior  domain  of  a  solid  continuum 


